//
//  MathHelper.h
//  SGD
//
//  Created by Niketan Pansare on 6/13/13.
//  Copyright (c) 2013 Rice University. All rights reserved.
//

#ifndef __SGD__MathHelper__
#define __SGD__MathHelper__

#include "Configuration.h"
#include <sstream>

// ===================================================
// Configuration parameters for user of this wrapper
// This parameters are not set in Configuration.h
// -------------------------------------------
// For sureshot answer (no reliability guarantees), Uncomment all these parameters:
// #define PRECISE_RATHER_THAN_FAST_OPTIMIZATION
// #define USE_MINUIT_OPTIMIZER
// #define USE_GSL_OPTIMIZER
// #define USE_LAST_ATTEMPT_OPTIMIZER
// -------------------------------------------

// Error handling macro: User can change this code depending on his application.
// #define HANDLE_GSL_ERROR(msg) { std::cout << "\nError: " << msg << "\n"; exit(1); }

// ===================================================

/// Installation step for Minuit and GSL (if you are using Ubuntu):
/// sudo apt-get install libgsl0-dev gsl-bin g++
/// wget http://www.cern.ch/mathlibs/sw/5_28_00/Minuit2/Minuit2-5.28.00.tar.gz
/// tar -xzf Minuit2-5.28.00.tar.gz
/// cd Minuit2-5.28.00
/// sudo ./configure
/// sudo make install
/// sudo ln -s /usr/local/lib/libMinuit2.so.0 /usr/lib/libMinuit2.so.0

// Compile:
// g++ -g myOtherCPPCode.cpp MathHelper.cpp `root-config --cflags` `root-config --libs` -lMinuit2 -fopenmp -lgsl -lgslcblas -lm -lstdc++ -o spoken_web

// Standard header files
#include <iostream>
#include <sstream>
#include <string>
#include <cmath>
#include <cassert>
#include <vector>
#include <cstring>


// GNU Scientific Libraries' header files
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_blas.h>

// Use Minuit header files if the flag is enabled
#ifdef USE_MINUIT_OPTIMIZER
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinimize.h"
#include "Minuit2/FCNBase.h"
#include "Minuit2/MnScan.h"
#endif

#define IS_FINITE(val) (gsl_finite(val) == 1)
#define LOG(val) ({if(val < 0) { std::cout << "Error: Incorrect parameter for LOG(" << (val) << ") at " << __FILE__ << ":" << __LINE__ << "\n"; exit(-1);} log(val);})
#define ARE_DOUBLE_EQUAL(val1, val2) (fabs((val1) - (val2)) < 1e-6)

// Note: SOFT_BOTH bounds are allowed to be wrong if you give good enough initial value.
enum TypeOfBoundsForOptimization { HARD_BOTH_SIDES, HARD_LOWER, HARD_UPPER, SOFT_BOTH };

#ifdef USE_MINUIT_OPTIMIZER
  class OneDimensionalOptimizer : public ROOT::Minuit2::FCNBase {
#else
  class OneDimensionalOptimizer {
#endif
  public:
    // ---------------------------------------------------------------------------------------------
    // Sample Usage:
    // std::vector<double> initialValues; initialValues.push_back(35);
    // OneDimensionalOptimizer myOpt(TheFunctionToBeOptimized, initialValues, 0.001, myLower, myUpper, SOFT_BOTH, 0);
    // double myMaxima, myArgMax;
    // if(!myOpt.FindMaxima(&myMaxima, &myArgMax)) {
    //  std::cout << "\nError: The wrapper class OneDimensionalOptimizer could not find a maxima\n";
    //  exit()
    // }
    // else {
    //  Yaay, maxima found .... USE myMaxima, myArgMax variables !!!
    // }
    
    OneDimensionalOptimizer(double(*myFunction)(const std::vector<double>& parameters), std::vector<double> initialValues, double tolerance, double lower, double upper, enum TypeOfBoundsForOptimization howBounded, unsigned int indexOfParameterToBeOptimized = 0);
    
    // Returns false if unable to find maxima of the function passed in the constructor
    bool FindMaxima(double* maxima, double* argmax);
    bool FindMinima(double* minima, double* argmin);
    
    // After FindMaxima, for rejection sampler only
    void FindTighterLRBounds(double maxima, double argMax, double* lower, double* upper, bool isLog);
    
    // Returns root: such that
    // 1. arg_low <= root <= arg_high
    // 2. rhsValue - tolerance <= f(root) <= rhsValue + tolerance
    double FindRoot(double rhsValue, double arg_low, double arg_high, double tolerance);
    
    // Binary search with tolerance. This still needs improvement. Ignore this for now.
    bool BinarySearch(double arg_low, double val_low, double arg_high, double val_high, double checkValue, double tolerance, double* retVal, unsigned int iter);
    
    // ---------------------------------------------------------------------------------------------
    
    // ---------------------------------------------------------------------------------------------
    // Ignore these, only implementor of this class has to worry about them.
    double operator()(const std::vector<double>& par) const;
    double Up() const;
    // ---------------------------------------------------------------------------------------------
  private:
    double(*f)(const std::vector<double>& par);
    bool doMinimization; double err;
    enum TypeOfBoundsForOptimization howBounded; double l; double r; unsigned int parameterIndex;
    bool doOptimization(double* value, double* argVal);
    std::vector<double> initialValues;
#ifdef USE_GSL_OPTIMIZER
    bool gsl_doesIntervalContainMaxima();
    bool gsl_oneDimensionalMaximization(double* value, double* argVal);
#endif
    double AbsoluteForDouble(double val);
  };
  
  // ---------------------------------------------------------------------------------------------
  // Other useful functions:
  // - Statistical functions: Mean (mean), Variance (variance)
  // - Log of special functions:  Gamma (lgamma), Factorial (lnfactorial)
  // - Absolute (AbsoluteForDouble)
  // ---------------------------------------------------------------------------------------------
  class StatisticsHelper {
  public:
    StatisticsHelper();
    ~StatisticsHelper();
    
    // --------------------------------------Begin Probability distribution -----------------------------------------------------
    /// \defgroup StatDistribution R-like API for family of probability distributions
    /// Naming convention:
    /// I have followed the naming conventions of R programming language. Please see its documentation for further details
    ///
    /// Here is a brief intro:
    /// Every distribution that R handles has four functions. There is a root name, for example, the root name for the normal distribution is norm. This root is prefixed by one of the letters
    ///   - d for "density", the density function
    ///   - r for "random", a random variable having the specified distribution
    /// Example: For the normal distribution, these functions are dnorm, and rnorm.
    ///
    /// Important note: Since most bayesian calculation is done in log space, (almost) every density function is capable of outputing 'numerically robust' logarithmic value. Hence, use dnorm(x, mean, sd, true) in your code rather than log(dnorm(x, mean, sd))
    ///
    /// Suported continuous distribution:
    /// - Uniform (unif - r,d)
    /// - Univariate Normal (norm - r, d); Truncated normal (tnorm - r,d);
    /// - Multivariate Normal (mnorm - r, d, conditional_r, conditional_rt) (--> Only m = 2, 3 supported)
    /// - Beta (beta - r,d)
    /// - Inverse Gamma (invgamma - r,d)
    ///
    /// Suported discrete distribution:
    /// - Discrete Uniform (dunif - r,d)
    /// - Dirichlet (dirichlet - r,d)
    /// - Multinomial (multinom - r,d)
    /// - Poisson (pois - r,d)
    /// @{
    
    /// @name Continuous Uniform Distribution
    ///@{
    /// Returns uniform variable [min, max)
    double runif(double min, double max);
    double dunif(double x, double min, double max, bool returnLog = false);
    ///@}
    
    /// @name Discrete Uniform Distribution
    ///@{
    /// Returns uniform variable [min, max) --> i.e. [min, max-1]
    int rdunif(int min, int max);
    double ddunif(int x, int min, int max, bool returnLog = false);
    ///@}
    
    /// @name Univariate Normal Distribution
    ///@{
    double rnorm(double mean, double sd);
    double dnorm(double x, double mean, double sd, bool returnLog = false);
    ///@}
    
    /// @name Univariate Truncated Normal Distribution
    ///@{
    double rtnorm(double mean, double sd, double lower);
    double dtnorm(double x, double mean, double sd, double lower, bool returnLog = false);
    ///@}
    
    /// @name Multivariate Normal Distribution
    ///@{
    /// Returns m numbers retVal[0] to retVal[m - 1]. The user is expected to allocate/deallocate the memory for these values in calling function.
    /// m - dimensionality (only 2 and 3 supported as of now)
    /// rhos[] = {rho} or rhos[] = {rho_xy, rho_xz, rho_yz}
    /// retVarStr can be "x" or "z", ...
    void rmnorm(double* retVal, double* means, double* sds, double* rhos, int m);
    double dmnorm(double* vals, double* means, double* sds, double* rhos, int m, bool returnLog = false);
    /// Conditional multivariate normal - All variables except the one specified in 'retVarStr' are given
    double cond_rmnorm(double* vals, double* means, double* sds, double* rhos, int m, const char* retVarStr);
    /// Conditional truncated multivariate normal (Similar to cond_rmnorm, except the returned value is sampled from truncated normal using the variable 'lower')
    double cond_rtmnorm(double* vals, double* means, double* sds, double* rhos, int m, const char* retVarStr, double lower);
    
    /// General multivariate normal (arbitrary m)
    /// mu is m x 1 vector and sigma is rowwise m x m matrix 
    void rmvnorm(double* retVal, double* mu, double* sigma, int m);
    double dmvnorm(double* x, double* mu, double* sigma, int m, bool returnLog = false);
    
    // TODO:
    // Something like this double cond_rtmvnorm(double* x, double* mu, double* sigma, int m, const char* retVarStr, double lower);
    ///@}
    
    /// @name Wishart and Inverse Wishart Distribution
    ///@{
    /// The user is expected to allocate/deallocate the memory for these values in calling function.
    /// dof > numDimension + 1 and scale is positive definite (so is the returned matrix).
    /// The matrices are in row-wise format.
    void riwishart(double* retVal, const unsigned int dof, const double* scale, int numDimension);
    void rwishart(double* retVal, const unsigned int dof, const double* scale, int numDimension);
    ///@}
    
    /// @name vonMises Fisher (vmF) distribution
    ///@{
    void rvmf(double* retVal, double* mean, double concentration, int numDimension);
    void rspherical_uniform(double* retVal, int numDimension);
    ///@}
    
    /// @name Beta distribution
    ///@{
    double rbeta(double shape1, double shape2);
    double dbeta(double x, double shape1, double shape2, bool returnLog = false);
    ///@}
    
    /// @name Inverse Gamma distribution
    ///@{
    double rinvgamma(double shape, double scale);
    double dinvgamma(double x, double shape, double scale, bool returnLog = false);
    ///@}
    
    
    /// @name Inverse Gamma distribution
    ///@{
    double rgamma(double shape, double scale);
    double dgamma(double x, double shape, double scale, bool returnLog = false);
    ///@}
    
    /// @name Dirichlet distribution
    /// Example: double alpha[] = {1, 1, 1}; double retVal[3]; rdirichlet(3, alpha, retVal);
    ///@{
    void rdirichlet(int sizeOfAlpha, double* alpha, double* retVal);
    double ddirichlet(int sizeOfAlpha, double* alpha, double* x, bool returnLog = false);
    ///@}
    
    /// @name Multinomial distribution
    /// Both prob and retVal should be of length sizeOfProb
    /// size, say N, specifying the total number of objects that are put into K (or sizeOfProb) boxes in the typical multinomial experiment.
    /// If the array prob is not normalized then its entries will be treated as weights and normalized appropriately. Furthermore, this function also allows you to provide log probabilities (in which case, you need to set isProbInLogSpace to true)
    /// Eg: double prob[3] = {0.1, 0.2, 0.7}; double retVal[3]; rmultinom(1, 3, prob, retVal, false);
    ///@{
    void rmultinom(int size, int sizeOfProb, double* prob, unsigned int* retVal, bool isProbInLogSpace = false);
    double dmultinom(int sizeOfProb, double* prob, unsigned int* x, bool isProbInLogSpace = false, bool returnLog = false);
    ///@}
    
    
    /// @name Poisson distribution
    ///@{
    int rpois(double lambda);
    double dpois(int x, double lambda, bool returnLog = false);
    ///@}
    
    // ----------------------------------------------------------------------------------------------------
    /// @name Rejection Sampler for parameters of special distributions: rParameter<TypeOfDistribution>(...)
    /// - Bivariate normal (2Norm - r), Trivariate normal (3Norm - r -> Sorry, this is not robust yet)
    ///
    /// Note:
    ///   - Depending on type of distribution parameterType can be: "mean_x", "mean_y", "mean_z", "sd_x", "sd_y", ...
    /// - x, y, z, mean_x, mean_y, ... are  arrays of length numValues
    ///   - For Trivariate, we assume Inverse gamma priors for mu's and sigma's and Beta prior for rho
    ///
    /// This case appears a lot in OLA paper and having this function as part of this class
    /// reduces lot of redundant code. For more details, contact Niketan.
    /// \n Also, in this case, we assume Inverse gamma priors for mu's and sigma's and Beta prior for rho
    ///@{
    double rParameter2Norm(double* mean_x, double* mean_y, double* sd_x, double* sd_y, double* rho, double* x, double* y, unsigned int numValues, const char* parameterType, double initialValue);
    double rParameter3Norm(double* mean_x, double* mean_y, double* mean_z, double* sd_x, double* sd_y, double* sd_z, double* rho_xy, double* rho_xz, double* rho_yz, double* x, double* y, double* z, unsigned int numValues, const char* parameterType, double initialValue);
    ///@}
    // ----------------------------------------------------------------------------------------------------
    
    ///@}
    // --------------------------------------End Probability distribution -----------------------------------------------------
    
    
    /// @name Special functions
    ///@{
    /// Returns arithmetic (or population) mean (i.e. (\sum_i array[i]) / n)
    double mean(double* array, int numOfElements);
    double variance(double* array, int numOfElements, double mean);
    
    double lgamma(double x);
    double lnfactorial(int val);
    ///@}
    
    /// @name Useful functions for operations in logscale
    ///@{
    /// Returns log(a' + b') where a = log(a') and b = log(b')
    double log_add(double a, double b);
    void NormalizeInLogScale(int sizeOfProb, double* logProb);
    ///@}
    
    double max(double a, double b);
    double min(double a, double b);
    
    double AbsoluteForDouble(double val);
    
  private:
    gsl_rng* global_random_number_generator;
    int AreRhosValidForTriVariateNormal(double rhoXY, double rhoYZ, double rhoXZ);
    int DoubleCompare(const void* a, const void* b);   // For sorting
    
    double rho_yz; double rho_xz;
    //int IsMatrixPositiveDefinite(gsl_matrix* A); -> Use matrix helper
    
    double GetMuBar(double mean_x, double mean_y, double mean_z, double sd_x, double sd_y, double sd_z, double rho_xy, double rho_xz, double rho_yz, double y, double z);
    
    double SampleUniformlyFromIntervals(double* intervalL, double* intervalR, int numIntervals);
    int GetIntervalsThatSatisfiesInequalityConditions(double* rhoXZs, double* rhoYZs, int numRhos, double* intervalL, double* intervalR);
    void updateIntersection(double* intervalL, double* intervalR, int numIntervals, double lower, double upper);
    void FindIntervalsForQuadraticInequalities(double a, double b, double c, double* lowerBound, double* upperBound, int* numIntervals);
    
  };
  
  class MatrixHelper {
  public:
    // --------------------------------------Begin Matrices -----------------------------------------------------
    /// \defgroup MatrixHelper Helper functions for matrices stored as 'row-wise' double array.
    /// I prefer this approach rather than create a seperate class.
    
    /// @{
    
    /// Returns determinant of matrix m of the given order
    double det(double* m, int order);
    
    /// Returns the inverse of matrix m of the given order. Remember the memory for retVal is allocated by the caller, not by this function.
    void inv(double* m, double* retVal, int order);
    
    /// Returns the inverse of matrix m of the given order. Remember the memory for retVal is allocated by the caller, not by this function.
    /// Note: Beware ... can you get away without using this ? See http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix
    void inv_LU_decomposition(const double *m, double *retVal, int order); // more general
    void inv_cholesky_decomposition(const double *m, double *retVal, int order); // works only for positive definite and is computationally faster than LU
    
    void outerProductOfVectorsToCreateMatrix(double* vector1, double* vector2, double* retMatrix, int numDimensionsOfVector);
    
    // Create matrix using: double *matrix = new double[numRows * numColumns];
    double get(const double* m, int i, int j, int columnWidth);
    void set(double* m, int i, int j, int columnWidth, double val);
    
    /// These functions expect the memory to be allocated by the caller
    void convert_rowWise_to_gslMatrix(const double* src, gsl_matrix * dest, int numRows, int numColumns);
    void convert_gslMatrix_to_rowWise(const gsl_matrix* src, double* dest, int numRows, int numColumns);
    
    // Return false if atleast 1 elements is not finite, else true
    bool areAllElementsFinite(const double* m, int numRows, int numColumns);
    bool areAllElementsFinite(const gsl_matrix* m, int numRows, int numColumns);
    
    bool isPositiveDefinite(const double* m, int numRows, int numColumns);
    bool isPositiveDefinite(const gsl_matrix* m, int numRows, int numColumns);
    
    void copy(const gsl_matrix* src, gsl_matrix* dest, int numRows, int numColumns);
    
    // Row-wise
    // Eg: cout << getString(m, 10, 10, " ", "\n");
    std::string getString(const gsl_matrix* m, int numRows, int numColumns, std::string columnSep, std::string rowSep);
    std::string getString(const double* m, int numRows, int numColumns, std::string columnSep, std::string rowSep);
    
    // zero, identity (1,1,1, ...), unit (0,0,1,0,0,...),
    // add, subtract, multiply, scalar_multiply (i.e. scale)
    // find, apply, equal,
    // transpose
    // eigen, diag,
    // is_null, is_pos, is_pos_definite, ...
    /// @}
    
    // --------------------------------------End Matrices -----------------------------------------------------
  };
    
  // Wanted to control the vector inside GradientDescent class, so created a helper class rather than a MyVector wrapper with similar overloaded functions.
  class VectorHelper {
  public:
    // Performs vector1 = vector1 + vector2
    void add(double* vector1, double* vector2, int sizeOfVector);
    double dot(double* vector1, double* vector2, int sizeOfVector);
    
    double angle(double* vector1, double* vector2, int sizeOfVector, bool in_degrees = false); // Find angle between two vectors
    
    void multiplyByConstant(double* vec, double constant, int sizeOfVector);
    
    // Performs retVec = vec1 + constant*vec2
    void multiplyByConstantAndThenAdd(double* retVec, double* vec1, double* vec2, double constant, int sizeOfVector);
    
    double absoluteValue(double* vec, int sizeOfVector);
    
    // vecTo <- vecFrom
    void copy(double* vecTo, double* vecFrom, int sizeOfVector);
    
    void createUnitVector(double* retUnitVector, double* vector1, int sizeOfVector);
    void convertToUnitVector(double* vec, int numDimension);
    
    // The caller has to call the delete[] method
    double* createEmptyVector(int sizeOfVector, double defaultValue = 0);
    
    void assignDefaultValue(double* vec, int sizeOfVector, double defaultValue = 0);
    
    void print(double* vec, int sizeOfVector);
    void printWithConstantMultiplier(double* vec, double constant, int sizeOfVector);
    std::string getPrintString(double* vec, int sizeOfVector);
    
    /// These functions expect the memory to be allocated by the caller
    void convert_rowWise_to_gslVector(const double* src, gsl_vector * dest, int sizeOfVector);
    void convert_gslVector_to_rowWise(const gsl_vector* src, double* dest, int sizeOfVector);
    
  private:
    // Returns log(a' + b') where a = log(a') and b = log(b')
    double log_add(double a, double b);
  };
    
  
  double NormalizeBetaToRho(double betaVal);
  double NormalizeRhoToBeta(double rhoVal);
  double myLogBetaPrior(double x);
  double myLogIGPrior(double x);

#endif /* defined(__SGD__MathHelper__) */
